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Abstract 

We present a brief review of the classical density functional theory of 
atomic and molecular fluids. We focus on the application of the theory 
to the determination of the solvation properties of arbitrary molecular 
solutes in arbitrary molecular solvent. This includes the prediction of the 
solvation free energies, as well as the characterization of the microscopic, 
three-dimensional solvent structure. 


1 Introduction 

The determination of the solvation free-energy of molecular solutes in molecular 
solvents is a problem of primary importance in physical chemistry and biology. 
From a theoretical point of view, two extreme strategies can be found in the lit¬ 
erature. A standard route consists in using molecular simulation techniques such 
as molecular dynamics (MD) or Monte-Carlo (MC), with an explicit molecular 
solvent. There are a number of well-established statistical mechanics techniques 
to estimate absolute or relative free-energies by molecular simulations, for ex¬ 
ample thermodynamic integration methods based on umbrella sampling [^|^, 
or generalized constraints . In any case, the precise estimation of free-energies 
by computer simulation remains extremely costly; it requires to consider a suf¬ 
ficiently large number of solvent molecules around the molecular solute and, 
for this large system, to average a "generalized force" over many microscopic 
solvent configurations, and this for a lot of different points along the reversible 
thermodynamic integration path. 

Another class of methods, known as implicit solvent models |^, relies on 
the assumption that the macroscopic laws remain valid at a microscopic level, 
and that solvation free energies can be computed by combining a dielectric con¬ 
tinuum description of the solvent outside the solute core and a simple solvent- 
accessible surface area expression for the non-electrostatic contributions (^. For 
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the electrostatic part, the stationary Poisson-Boltzmann equation can be solved 
for the electrostatic potential using sharp definitions of the dielectric bound¬ 
aries and various efficient numerical techniques, making it possible to handle 
very large biomolecular systems [^. Density functional methods based on the 
minimization of polarization density have been introduced too. There are 
serious limitations however to a continuum dielectric approach, and first of all 
the validity of the macroscopic electrostatic laws at microscopic distances, the 
neglect of the molecular nature of the solvent, and the ambiguous definition of 
all the non-electrostatic energetic contributions, such as hydrophobicity. Macro¬ 
scopic approaches to hydrophobicity that mixes consistently with the Poisson- 
Boltzmann description are presently developed. 

Beyond continuum descriptions, it is desirable however to devise and employ 
implicit solvent methods which (i) are able to cope with the molecular nature 
of the solvent, but without considering explicitly all its instantaneous micro¬ 
scopic degrees of freedom, and (ii) can provide solvation properties at a modest 
computer cost compared to explicit simulations. Such methods should rely on 
the theory of molecular liquids that has been developed in the second half of 
the last century and lie by now in classical textbooks 10 -12 . Among possible 


16 or molecular 17-23 , or mixed 1241^ 


approaches one should mention molecular integral equation theories in the ref¬ 
erence interaction site (RISM) 
picture, Gaussian field theories 
molecular liquids [Toj 28 


2^ 27 , the density functional theory (DPT) of 


or. 


finally, field theoretical approaches to dipolar 
solvent-ions mixtures, that lead to a generalization of the Poisson-Boltzmann 
equation accounting for particle size and dielectric saturation j^-44 . Note that. 


close to our purpose, a 3D-version of the RISM equations has been developed 
recently to describe the solvation of objects of complex shape [45 -50 


Our focus here is classical density functional theory (DPT), and eventually 
a molecular version of it that we will call Molecular Density Punctional The¬ 
ory (MDPT). The basic theoretical principles of classical DPT can be found 
in the seminal paper by B. Evans |28| and subsequent excellent reviews by 
him (28}j30] and other authors 51 . The advent in the late 1980’s of a quasi- 


exact DPT for inhomogeneous hard sphere mixtures, the fundamental measure 
theory [^ - 57 , has promoted recently a great deal of applications to atomic- 
like fluids in bulk or confined conditions or at interfaces. Classical "atomic" 
DPT can be considered nowadays as a method of choice for many chemical en¬ 
gineering problems 58,59 . Much less applications exist for molecular fluids, 
for which solvent orientations should be considered. The description has been 
generally limited to generic dipolar solvents 60 61 or dipolar solvent/ions mix¬ 
tures 33 -36 ; such approach may be considered already as "civilized" compared 
to primitive continuum models |35j . We have proposed recently an extension of 
MDPT to arbitrary fiuid/solvents in the precise goal of describing the solvation 
of three-dimensional molecular object in arbitrary solvents. (37-39,62-^ A 


RISM-based DPT approach of molecular solvation has been developed recently 
too (tI) . 

The outline of the present review is as follows. We first recall the basic 
principles of cDPT for atomic-like fluids. We then describe the particular but 
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fundamental case of the hard-sphere fluid, and the associated fundamental mea¬ 
sure theory (FMT), focusing on a "scalar" formulation due to Kierlik and Ros- 
inberg (^[5^ , instead of the standard "vectorial" version introduced initially 
by Rosenfeld [^. We then turn to Lennard-Jones fluids, for which the HS 
fluid can be used as a reference in various ways to construct a functional. The 
last section will be devoted to molecular solvent, modeled by rigid polyatomic 
molecules with an orientation. The discussion will focus on a model dipolar sol¬ 
vent, the Stockmayer fluid, and then extend to realistic models of polar liquids 
such as acetonitrile and water. 


2 The case of atomic fluids 


2.1 General formulation 

In this section we begin with recalling the basis of the density functional theory 
of liquids, and discussing the general problem of a molecular solvent submitted 
to an external held. In the applications we have in mind, the external held will 
be created by a molecular solute of arbitrary shape dissolved at infinite dilution 
in the solvent. The individual solvent molecules will be considered later as rigid 
bodies described by their position r and orientation oj. In this section we restrict 
the discussion to atomic or pseudo-atomic solvents (such as CCI 4 ) modeled by 
spherical particles for which only the position r matters. 

The grand potential density functional for a fluid having an inhomogeneous 
density p{v) in the presence of an external held V^xt (r) can be defined as 28p9 


^[P] = F[P\ ~ j 


( 1 ) 


where F[p\ is the Helmholtz free energy functional and ps is the chemical poten¬ 
tial. The grand potential can be evaluated relatively to a reference homogeneous 
fluid having the same chemical potential /is and particle density po 


H[p] = Vt[pQ\+F[p]. 


( 2 ) 


Following the general theoretical scheme introduced by Evans 10 28 29,72 


the density functional F[p\ can be split into three contributions: an ideal term, 
an external potential term and an excess free-energy term accounting for the 
intrinsic interactions within the fluid, 


T F" [p] -j- J~ ea:c[p]5 

with the following expressions of the first two terms 

Xri) 


p(ri)ln 


J^id[p\ = ksT J dr 1 
Fext[p] = J driVext{ri)p{ri). 


Po 


-p(ri) -kpo 


( 3 ) 

( 4 ) 

( 5 ) 
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There are several ways to arrive to an exact expression of the excess free-energy, 
i.e. using an adiabatic perturbation of the pair potential (the so-called adiabatic 
connection route in electronic DFT), of the external potential, or of the density 
itself. If the latest route is chosen, one can define J-exc[p] ^-s: 

^exc[p\ = ksT JJ dridr2C{ri,r2) Ap{ri)Ap{r2), ( 6 ) 


with Ap (r) = p{r) — po. The function C'(ri,r 2 ) is still a functional of p (r) 
defined by 


C'(ri,r 2 )= [ da{a-l) c^^^(ri, r 2 ; [pa]), 
Jo 


(7) 


where c^^^(ri, r 2 ; [pa]) is the two particle direct correlation function, i.e. by 
definition the second order derivative of the excess free-energy with respect to 
density, evaluated at the intermediate density pQ(r) = po + aAp(r). Eqs 
follows naturally when expressing the functional from the knowledge of its 
second-derivatives 


28 


The equilibrium condition reads 

5n[p] 


Sp 


= 0 


5J-[p] 


Sp 


= 0 . 


( 8 ) 


When minimizing the density functional JF[p] with respect to p (r), the value at 
the minimum is the difference of the solvent grand potential with and without 
the solute, and thus the solute solvation free-energy. The associated density 
Pe 5 (r) is the equilibrium inhomogeneous density. 

The functional defined by eqs [3][^ is formally exact but the inhomogeneous 
direct correlation functions entering the definition of the excess term are indeed 
unknown. However, simple approximations can be proposed for this quantity. 
The most natural one consists in expanding the inhomogeneous direct correla¬ 
tion function c*'^^(ri, r 2 ; [pa]) around a = 0, that is, around the homogeneous 
density po: 


(r 1 , r 2 ; [pa]) = (ri, r 2 ; [po]) -f a 


dr 3 c(^)(ri,r 2 ,r 3 ; [po]) Ap(r3)-k.... (9) 


Such expression involves the two, three, n-particle direct correlation functions of 
the homogeneous fluid [^. The first term is the (two-body) direct correlation 
function (DCF) of the homogeneous solvent, that depends on ri 2 = |r 2 — ri|, 
and can be thus denoted as C 5 (ri 2 ;po) {S for spherical component, preparing 
ourselves to non-spherical solvents). 

Using eq. the excess term can thus be written as 

J^exc[p\ = -^ dridr2Cs{ri2]Po) Ap{ri)Ap{r2) + J^bIp], ( 10 ) 

where we have defined the bridge functional 

J^b[p] = ~^ j dridr2dr3 c^^\ri,r2,r3;[po]) Ap{ri) Ap{r2) Ap{r3) + o{Ap*), 

( 11 ) 
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which thus starts with a cubic term in Ap. Setting Tb[p\ = 0 correspond to the 
so-called homogeneous reference fluid (HRF) approximation. It can be shown 
to be equivalent to the hypernetted chain (HNC) approximation in integral 
equation theories 30 . The input of the theory is the direct correlation function 


of the pure solvent, which can be extracted from simulation or experimental 
data by measuring the total correlation function hs(r) = g{r) — 1 and solving 
subsequently the Ornstein-Zernike equation, i.e. in Fourier space: 


I-Pocsik) = {1 +pohsik)) =x„ (k). 


( 12 ) 


Xni'i') is the structure factor, or the density susceptibility, measuring density- 
density correlations at a given distance in the fluid. The excess free energy can 
thus be expressed also in terms of the susceptibility 




2.2 


= [p] = J *idr2X„^(ri2) Ap(ri)Ap(r2) - J dr Ap{r f + J^b[p]- 

(13) 

Fundamental Measure Theory for the hard-sphere fluid: 
Scalar versus Vectorial Formulation 


We focus here on the particular case of the hard-sphere fluid and describe briefly 


the fundamental measure theory (FMT) introduced by Rosenfeld 52 and Kier- 
lik and Rosinberg j^. Although the theory is valid for arbitrary hard-sphere 
mixtures, we consider a one-component HS fluid composed of hard spheres of 
radius R and at a bulk density po- The fluid is subjected to an external pertur¬ 
bation, for example a solid interface or a molecular solute of arbitrary shape em¬ 
bedded in the fluid, that creates a position-dependent external potential 14a:t(r) 
and thus an inhomogeneous density p{r). The excess functional of eq. [^can be 
written as 




Ap] = Flic [p] - Fiif [po] - PeJi J dr (p(r) - po ), ( 14 ) 

where (p) is the excess Helmholtz free-energy functional for the hard-sphere 
fluid and is the bulk excess chemical potential defined by 


..HS _ dF^^flp] 

r^exc 


dp 


P=Po ) 


(15) 


so that obviously 


ea;c[p] I 


o=po = 0. In the FMT introduced by Rosenfeld 


52 


Sp 'P — P^ 

the excess functional for the hard-sphere fluid can be written in terms of a set 
of Au, weighted densities, na(r): 


with 


Fexcl{pi(r)}J = kBT J dr$({na(r)}) 
Hair) = J dr'p{r')u}a{r-r') = p{r)-ku}a{r), 


(16) 

(17) 
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where (r) are geometrical weight functions to be defined below and ★ indicates 
the convolution of the microscopic densities by those weight functions. 

In the original Rosenfeld’s derivation there are four scalar weight functions, 
a;^(r), a = 0,1, 2,3, and two vectorial ones ci;i(r), a; 2 (r) that are defined by 


W3(r) 

II 

© 

1 

(18) 

^ 2 ( 1 ') 

= 47ri?wi(r) = 47ri?^ a;o(i’) = (5(R — 7-) 

(19) 

W2(r) 

= 47ri?d;i(r) = - ^(i? — r). 

(20) 


0(r) denotes the Heaviside function and 6{r) the Dirac distribution. The ex¬ 
cess free-energy density $ derived by Rosenfeld for Eq. [T^ is a function of 
the four position-dependent weighted densities, na{v),a = 0,1,2,3, and of the 
two vectorial ones, ni(r),n 2 (r), which generates in the homogeneous limit the 
Percus-Yevick equation of state for hard-sphere mixtures. Starting from the gen¬ 
eralization of the Carnahan-Starling (CS) equation of state to mixtures (namely 
the Mansoori-Carnahan-Starling-Leland equation (MCSL)) instead of PY, Roth 
and Wu et al 


et al 55 


56 were later able to obtain a modified expression based 


on the same definition of the weighted densities (either called white-bear (WB) 
version or modified FMT version (MFMT)). This modified version of FMT takes 
advantage of the fact that the CS expression provides one a better equation of 
state that PY. 

Ten years before those latest developments, Kierlick and Rosinberg were able 
to derive an alternative version of FMT which involves only four scalar weight 
functions a;a(r),a = 0,1,2,3. 54 . The last two weights are identical to 


Eq. [I8p9| whereas the first two ones are given by 


wi(r) = —6'{R-r) 

^o(r) = + 


( 21 ) 

( 22 ) 


Those weight functions appear naturally in the derivation as the inverse Fourier 
transforms of 


UJ3(k) 

Ul2(k) 

Wl(fc) 

Wo(fc) 


47r 

—risinikR) — kRcos(kR)) 
k^ 

AttR 

—-— sm{kRi) 


— (sm{kR) + kRcos{k)) 
2k 

kR 

cos{kR) + — sm{kR). 


(23) 


Although the main part of the papers by Kierlik and Rosinberg relies on a PY 
expression for the excess free energy density 


= -no ln(l - na) -f 

i — 713 


-k 


1 nl 

247 r (1 — 713)2 ’ 


(24) 
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the authors do mention in their conclusion that a CS (more precisely MCSL) 
expression could be used instead 




/J_n| 
\ SGtt n§ 


no) ln(l - no) + 

J 1 - no 


1 nl 

367r (1 - no)2no' 


(25) 


They point out the fact that this expression is more precise than the PY one, 
but using it while keeping the expression of the weights unchanged leads to 
thermodynamic inconsistencies; those inconsistencies are indeed present in the 
WB or MFMT formulations too. There is clearly a trade off to be made between 
precision and theoretical consistency. It was later shown by Phan et al. that the 
Kierlik and Rosinberg’s approach is mathematically equivalent to the original 
vectorial version. |73| On a practical point of view, however, and especially 
in the perspective of 3D applications, the KR formulation is advantageous with 
respect to the Rosenfeld’s formulation since the number of independent weighted 
densities is reduced from 5 to 4 for the one component system, and from (10 + 
Ng) to {4 + Ng) for a mixture of Ng components, with Ng > 2; thus from 
12 to 6 for a binary mixture. An efficient numerical implementation in three- 
dimensions of the Kierlik-Rosinberg FMT functional is detailed in Ref. [^. The 
numerical efficiency of the algorithms, in terms of convergence rate and system 
size dependency, is briefly illustrated in Fig. ^ 


2.3 The Lennard-Jones fluid 


Building the thermodynamics of the Lennard-Jones fluid by taking the suitable 
hard-sphere fluid as a reference and building in the attractive interaction as a 
perturbation is indeed a classics in liquid state theory and is at the basis of the 
Van der Waals theory of fluids. When coming to functionals, this idea can be 
declined in several variants, the most natural one being to use a FMT func¬ 
tional for the repulsive part and a mean-field approximation (or mean spherical 
approximation, MSA) for the attractive part [^. Along the lines given above, 
another route is to approximate the bridge functional in eq. [^by a hard sphere 
bridge functional, introduced by Rosenfeld as a universal bridge function |74l[75] 


^■[^(r)] = ksT J dr 


p{r )In 


Po 


- p{^) + Po 


-k J dr Vext{r)p{r) 


ksT 

2 


J dridr 2 cs{ri 2 \po)Ap{ri) Ap{r 2 ) + J^b^[p], 


(26) 


where 


= F^J[p{r)]-F^J[po]- J drAp{r) 

+ J dridr2cf®(ri2;po)Ap(ri) Ap(r2). (27) 


The first three terms represent the one-component hard-sphere KR-FMT excess 
functional defined in the previous section and the associated chemical potential 
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Figure 1: Top: Typical plot of the free energy difference between two successive 
steps (normalized by the initial energy) versus minimization-step number (Here 
a benzene molecule in a one-component HS reference fluid modeling SPC water). 
The inlet represents the same with a logarithmic scale in ordinates. Bottom: 
CPU time per minimization step versus number of 3D-grid points. The circle 
correspond in increasing order to N= 32, 64, 128, and 256. 










yielding equilibrium at p(r) = pQ. The fourth term involves the direct correlation 
function of the HS fluid at the same density, i.e 


ffS/i I \ M I 

cs (|ri r 2 |,po)- 

This function can be easily obtained in Fourier space as 


„HS 


ik;po) = , ({^ 7 }) u}a{k)ujp{k), 




drLadnp 


(28) 


(29) 


where {u^} represent the weighted densities for a uniform fluid of density po 


and the uja,p{k) are the weights of eq. 23 The second derivatives have to be 


taken for the PY or CS functions of eqs|24|or[2^ Note that defined as in eq. 
■^B^[p(^)] carries an expansion in Ap of order 3 and higher that corrects the 
second order expansion of the excess free energy in eq. The excess functional 
can also be re-expressed as 


= F^J[p]-FfJ[po]- J drAp{r) 


ksT 


dridT2 cf\ri2; Po) Ap(ri)Ap(r2), 


(30) 


where we have defined the "attractive" DCF by 

cT{ri 2 -,Po) = cs{ri 2 -:Po) - cf‘^(?’i 2 ; Po)- 


(31) 


Eq. 30 is the basis of the first order mean-spherical approximation (FMSA) 


theory developed by Tang 76 


We show here how this FMSA theory works for our purpose: the prediction 
of solvation properties of dissolved molecular objects. In Fig. |2.3| we compare 
the solvation free energy of a LJ sphere of increasing diameter in a LJ fluid with 
p* = 0.85, T* = 0.88, as computed by Monte-Carlo simulations by Lazaridis ff7\ , 
to the results obtained by DFT minimization with different HS diameters, d. 
It can be seen that the results are extremely sensitive to the choice of d, and 
that the best agreement is obtained for d = 1 .Older (indeed close to 1 , that 
would be the initial guess value). For that value, we have plotted in Fig. 2.3 
the microscopic solvent density, g{r) = p{r)/po, obtained for solute of different 
sizes by direct simulation, or by DFT in the HNC or FMSA approximation. It 
can be seen that the addition of hard-sphere bridge in FMSA greatly improve 
the results compared to the HRF (or HNC) approximation and yields a correct 
structure. 
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Figure 2: Solvation free-energy obtained by DFT using the HS bridge func¬ 
tional of with different HS diameters, compared to the Monte-Carlo resnlts 
of Lazaridis . 


3 The case of molecular fluids: Molecular deusity 
fuuctioual theory (MDFT) 

3.1 General formulation 


The solvent molecnles now carry a molecnlar structure that is described by a 
collection of distribnted atomic interaction sites. The theory is formulated in 
the molecnlar pictnre in which each solvent molecnle is considered as a rigid 
body and characterized by its position, r (e.g. the position of center of mass), 
and by its orientation, w, defined by the three Euler angles cu = (0, (j), ip). Thns, 
in the presence of an external perturbation, the solvent is now characterized by 
an inhomogeneous position and orientation density p{r,u}). The solute, as the 
solvent, is described in microscopic details by a molecular non-polarizable “force- 
field" involving atomic Lennard-Jones and partial charges parameters. Given 
that the solnte is fixed and defined by the position, Rj, of its different atomic 
sites, the external potential is defined by 


Fea:t(r,w) — ^ ^ 

i^solute j^solvent 




Aneorij ’ 


(32) 


where eij and aij are the Lennard-Jones parameters between solute site i and 
solvent site j, and qi and qj are the partial charges carried by those sites. 
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Figure 3: Reduced solvent density around LJ solutes of different diameters, 
using the HNC approximation, or adding a hard-sphere bridge functional with 
d = l.OMcr. 
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The relative site-site vectors are function of the solvent molecule position and 
orientation and defined as r^- = r -f R(a;)sj — R^, where Sj denotes the site 
positions in the molecular frame and R(a;) is the rotation matrix associated to 
10 . 

The same density functional as in eqs[l|^can be written for p(r,a;), with 
an ideal, external, and excess part: 


^id [p] 
^ext [p] 
^ex c [p] 


ksT / drdio 


p(r,w))ln 


Stt^P (r, io) 


no 


- p(r,uj)) + 


no 

Stt^ 


dTdujVext{r,uj)p{r,uj) 


(33) 


1. 


-2^bTJJ dridr2duJiduj2Ap{ri,ioi)c{ri-r2,uJi,i02)Ap{r2,uj2) 
+TB[p{r,uj)], 


where Ap (r, lo) = p (r, uj) — no /Stt^, uq being the particle number density of the 
reference bulk fluid. The first term represents the Homogeneous Reference Fluid 
approximation (or HNC approximation) where the excess free-energy density is 
written in terms of the angular-dependent direct correlation of the pure solvent. 
The second term represents the unknown correction to that term (or Bridge 
functional) that, again, can be expressed as of a systematic expansion of the 
solvent correlations in terms of the three-body, .. n-body terms direct correlation 
functions. 


3.2 The Stockmayer solvent 

To test and illustrate the theory, we start from the simplest conceivable model 
of dipolar solvent, the Stockmayer model, characterized by a single Lennard- 
Jones center with parameters crs,es ^^nd a dipole = pto, where lo is the 
unitary orientational vector of the molecule -which here replaces the orientation 
noted lo above. The parameters are selected to make the model look like water 
(similar density, no = 0.033 particles/^, particle size, Us = 3, and molecular 
dipole, p = 1.851?) although not tasting quite as water (no hydrogen bond in 
the model!). For such model, the external potential can be written as 


with 


Vext(r,(o) = $Lj(r) - pE,(r) • lo 


M 


$L.7(r) = 


Z=1 


|r - R, 


|r - R, 


M 


^ 4^eo^ |r-R,|3 ■ 


(34) 


(35) 

(36) 


It is also argued in Refs 19 that the c-functions can be expanded onto a 
rotational invariants basis set keeping, to a good approximation approximation. 
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the same order as the intermolecular potential, namely 

c(ri2,u;i,W2) = cs(ri2) +CA(ri2)«'“°(wi,W2)+ CD(ri2)$^^^(u;i,a;2037) 

where 


$ 

$1 


no 


= uji-u:2, 

= 3 (a;i • f 12 ) (u;2 • f 12 ) - wi • 0^2 


(38) 


represents the two first non-isotropic spherical invariants. The three compo¬ 
nents C 5 ,ca,C£i of c(ri 2 ,tiJi,(^ 2 ) can be obtained from the the corresponding 
components of the total correlation function, /i(ri 2 , Wi, 0 ^ 2 ), by inversion of 
the angular-dependent OZ equation. The total correlation function itself can 
be measured by using, e.g., MD simulations. We have used here the origi¬ 
nal Wertheim’s notations with subscripts S, A, and D for the different h- and 
c-components. In this approximation, it was shown in Refs 37 39 that the 


OZ equation can be solved directly for the different components in real space. 
Results of equivalent precision can be reached from inversion relations in k- 
space [TS] . 

Defining the number density. 


i(r) = J du}p{r,u}), 


and the polarization density. 


P(r) =/r J duj u) p{r, u)), 


(39) 


(40) 


and injecting the expression 37 of c(ri 2 ,u)i,u> 2 ) into the functional of eq. 33 it 


can be shown that external and excess terms can be written as functionals of 
n(r) and P(r) instead of the much more complex variable p{r,u}), namely 


J^extln^V] = 


exc[n,P] — 


y dr$Lj(r) - y drP(r) • Eg(r) 

y dvxdv2 cs(ri 2 )An(ri) • An(ri) - J dridr2 ca(?’ 12 )P(i'i) • P(ri) 


2 

keT 
' 2^2 


y dridr 2 cd(?' 12 ) [3(P(ri) • fi 2 ) (P(r 2 ) • fi 2 ) - P(ri) • P(ri)] . 


(41) 


At this stage, the expression of the ideal term can be kept unchanged as a 
function of p (r, w) as in eq. 33 and the minimization of the whole functional still 


performed with respect to p (r, w). The above expressions of the nonlocal excess 
free energy requires to perform FFT’s for n(r), P(r) rather than for p (r, w) and 
this reduces considerably the computation time. We can go even a little bit 
further, and show that the ideal part itself can be expressed as a functional of 
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n(r) and P(r) [^, namely 


^id [^1 P] 


r Tzfrl 

ksT / drn(r)ln(-) — n(r) + ng 

J no 

ksT J drn{r) I In 


(42) 


r-ltJHlL) 

^ yfin{r)) 




-P(r) P(t) 


^n(r) /in(r) 


In the second, polarization term, C designates the Langevin function and its 
inverse; P{r) is the modulus of the polarization vector P(r). The linearization 
of this term for small polarization fields yields the correct electrostatic limit, 
namely 

where ad = fj?/iksT is the usual equivalent polarizability of a dipole fi at the 
temperature T. One recognizes the expression of the polarization free-energy in 
a medium with local electric susceptibility Xe(r) = adn(r). 

Although the functional is now complete and usable as such, we proceed by 
looking at an equivalent of eq. involving susceptibilities rather than direct 
correlation functions. We introduce the longitudinal and transverse polarization 
in k-space 


PL(k) = (P(k).k)k 

Pr(k) = P(k)-Pi(k), (44) 


where k = k/fc. The electrostatic part of the excess free energy in eq 
written in k-space 


41 


can be 




■elec 

exc 


1 ksT 

2 ^2 

1 ksT I 

2 J 


J dkcA(A:)P(k) •P(-k) 

dkcoik) f3(P(k) • k)(P(-k) • k) - P(k) • P(-k)l (45) 


This can be easily rearranged into 


jrelec ^ 1 ksT 

exc 2 /i2 


dkc_(fc) PT(k)-PT(-k)-|- 


with the usual definitions llO 79,80 


dkc+(A:)PL(k).Pi(-k) 

(46) 


c-{k) = CA{k) - coik) (47) 

c+{k) = CA(fc) + 2cd(/c). (48) 


We use now the relations between C-{k) and c+(k) and the longitudinal and 
transverse dielectric constant eL{k) and eT^k), or, alternatively, the longitudinal 
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and transverse dielectric susceptibilities Xi(fc) and XTik) (see Refs [W 79-^) 

l-l/eL{k) 47rxL(fc)’ 

3y 3y 


l-yC+(A:) = 


l-^c_(A:) = 

with y = fi^no/9kBTeo, such that 


eT(fc) - 1 47rxT(fc)’ 


(49) 

(50) 


F. 


elec 


+ 


SksT 
noiF . 

1 

Sttco 


J dkP(k) • P(-k) 

■^kPT(k).P.(-k)^ 


dk 


Pi(k).PL(-k) 


(51) 


Xrik) J XL{k) 

(beware of the definition of xl, with or without a 47r factor |79||81|). 

At the end, one can gather all the above equations, including eqs. mm 
to get the following functional for a dipolar fluid, defined in terms of the density 
and dielectric susceptibilities 


Fid[n,'P] = ksT / drduj 


p (r, Lj) In 


47rp (r, Lj) 


riQ 




+ / dr$L, 7 (r)n(r) - 


fcaT 
2nn , 


dv An(r)^ 


keT 


dvidv2Xn (ri2) Ap(ri)Ap(r2) 
3k bT 


- drP(r).E,(r)- 


no/i^ 


drP(r)" 


(52) 


1 


Stteo 

1 

STren 


J dridr2 Xt^(^i2) Pt)!"!) 
J dridr2 xl\ri2)P Liri) 


PT(r2) 

PL(r2) +FB[n,P]. 


The ideal part can also be taken as in eq. so that the whole functional can 
be minimized with respect to n(r) and P(r). The bridge term can be neglected 
(HNC approximation) or approximated as a functional of n(r) only, FB[n\, for 
example using the hard-sphere bridge functional of Section |2.2| 

As a short illustration, Fig. Hr shows the accuracy of the MDFT approach 
(within the HNC approximation) for the microscopic structure of the Stock- 
mayer solvent around neutral and charged spherical solutes 38 63 . The 
MDFT results are compared to direct MD simulations of the solute embed¬ 
ded in the solvent. They do appear very satisfactory and account accurately 
for the shape of the peaks and their variation with charge and size (despite a 
slight overestimation of the first peak height for the neutral solute). Fig. & 
illustrates the case of a multisite polar molecule (here a three-site model of the 
acetonitrile molecule) with similar conclusions. An application to a more com¬ 
plex molecular system, namely the three-dimensional solvation structure close 
to an atomistically resolved clay, is illustrated in Fig. and described further 
in Ref. (^. 
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Figure 4: Top: Reduced density of the Stockmayer solvent around various so¬ 
lutes. MDFT results (solid black lines) are compared to MD simulation results 
(dashed red lines). From left to right: CH 4 , Cl“, K+. Bottom: Same than for 
the various sites of an acetonitrile molecule dissolved in the Stockmayer solvent. 
From left to right: CH 3 , C, N. 
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Figure 5: Two-dimensional maps of the solvent number density n(r)/no in three 
different planes close to a neutral clay surface, as calculated by molecular dy¬ 
namics (top) and HRF-MDFT (bottom). Those planes correspond to a prepeak 
(left), the first maximum (center) and second maximum (right) of the out-of 
plane mean solvent density. See Ref. 
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Figure 6: Ion-oxygen pair distribution functions for chloride and potassium 
ions in SPC/E water computed by MD (red lines) or MDFT without (blue 
lines) or with the three body term described in Ref. (black lines) 
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Figure 7: A selection of solute site-oxygen pair distribution functions for the 
n-methyl-acetamide molecnle CH 3 NHCOCH 3 (shown on top) computed by MD 
(red lines) or MDFT (black lines), including the three-body term described in 
Ref. 1^ . 
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3.3 Extension to water and arbitrary molecular solvents 


Water is a special case, certainly by its very subtle physics, but also for the 
fact that the most popular molecular models fall in the category of simple point 
charge models with a single Lennard-Jones center (usually centered on the oxy¬ 
gen atom) and distributed point charges. In that case, it was shown recently 
that the functional just displayed in eq. (with the linear orientation vector 
u> substituted by the three-angle orientation w, and 47r by Stt^) is perfectly ap¬ 
plicable if the dipolar polarization P(r) is replaced by a multipolar polarization 
vector, accounting for the full charge distribution of the water molecule, and 
defined in k-space by 

P(k) = J da;/i(k, a;)p(k, w) (53) 

with 


/x(k,w) 


E Sm(t^) 

^”"k- s^(a;) 


A k-i 



/x(w) -k 


^ ^ g™ (k • 


Sm(w)) Sm(w) + 


(54) 

(55) 


being defined as the polarization, of a single molecule located at the origin. 
Sm(fl) designates the location of the atomic site for a given orientation 
uj. It reduces to the usual molecular dipole /x(a;) = 9m Sm(w) at dominant 
order in k. The multipolar dielectric susceptibilities xl{X) and xtA) entering in 
eq.|^can be either computed from MD simulations of the pure liquid, according 
to the procedure in Refs 82 , or inferred from experiments. 

Although giving already sensible results for rather complex systems [^, it 
was shown in that the HNC approximation Tb = 0 turns unfortunately short 
for describing the solvation of hydrophobes , as well as that of molecular 

solutes giving raise to strong H-bonds . Three-body corrections, including a 
spherical HS bridge, or a three-body term re-enforcing tetrahedral order, have 
to be added to give correct solvation structure and thermodynamics. This is 
illustrated in Fig. for the hydration structure around monovalent ions. Fig. 
shows the water structure obtained by MDFT around a N-methyl-acetamide 
molecule (the prototype for a NH-CO peptide motif), including the three-body 
correction term in the functional. |68| 

For a general solvent with more complex geometry, and described by more 
than one Lennard-Jones center, the full angular-dependent functional of Sec. |3.1| 
has to be adopted, and the necessary input remains the full angular-dependent 
direct correlation function c(ri 2 ,Wi,u; 2 ). Remaining in the HNC approxima¬ 
tion, this formulation was applied with some success to the study of charge 
transfer processes in acetonitrile (^ . 

The MDFT approach is still under current development, as are related site- 
DFT approaches [^, for practical applications such as the systematic prediction 
of solvation free energies . Classical density functional theories are expected 
to provide soon an alternative to the 3D-RISM approach, which is nowadays 
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becoming quite popular for applications in biological and material sciences - 
despite some intrinsic theoretical limitations that specialists are aware of. 
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